We’ll carry on looking at spatial data in this lab. You should have the following files from the previous lab, but it’s worth checking to make sure you have everything:

  • Climate dataset for Western North America: WNAclimate.csv
  • Temperature dataset for Oregon in shapefile format: oregon.zip
  • New York state polygon data in a shapefile: NY_Data.zip
  • Digital elevation map from Switzerland: swiss_dem.grd
  • A NetCDF file of global monthly air temperature: air.mon.ltm.nc
  • A NetCDF file of global monthly mean air temperature: air.mon.mean.nc
  • A dataset of country socioeconomic variables: countries.zip

Base R has no structure for spatial data, so you will need to install the following packages (you should have some of these from previous modules):

  • sf
  • terra
  • RColorBrewer
  • ggplot2
  • viridis
library(ggplot2)
library(terra)
library(RColorBrewer)
library(sf)
library(viridis)

Raster data

Previously, we were working with vector spatial data. These are geometries composed of points defined by their coordinates. An alternative form of spatial data is known as a raster. This is gridded data. It takes the form of a rectangle composed of squares of equal size, which are sometimes called ‘cells’ or ‘pixels’. Each cell stores some kind of value.

This simplifies the geometry, which can be specified by two pieces of information: the spatial extent of the raster and the resolution of the cells. Here we create a blank raster with 10 rows and 10 columns, with a resolution of 10x10 using the rast() function. We then assign random values to each cell:

r <- rast(nrow = 10, ncol = 10, 
          xmin = 0, xmax = 100, 
          ymin = 0, ymax = 100)

r[] <- runif(n = 100)

r
## class       : SpatRaster 
## dimensions  : 10, 10, 1  (nrow, ncol, nlyr)
## resolution  : 10, 10  (x, y)
## extent      : 0, 100, 0, 100  (xmin, xmax, ymin, ymax)
## coord. ref. :  
## source(s)   : memory
## name        :       lyr.1 
## min value   : 0.004815464 
## max value   : 0.972897647

rast objects can be plotted using the base plot() command:

plot(r)

The terra package offers a wide array of functions for dealing with gridded data, including the ability to read from many widely used file formats, like remote sensing images (e.g. GeoTiffs), NetCDF, and HDF formats. We will use it here to work with a Landsat 8 scene collected on June 14, 2017. The subset covers the area between Concord and Stockton, in California, USA. The files are contained in the zip file rs.zip. Download this (if you haven’t already) and unzip it in your data folder.

We will also need the shapefile of California places (ca_places.zip), so download and unzip this as well. We’ll read this in before starting:

ca_places <- st_read("./data/ca_places/ca_places.shp")
## Reading layer `ca_places' from data source 
##   `/Users/u0784726/Dropbox/Data/devtools/ugic-intro-r/data/ca_places/ca_places.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1618 features and 16 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -124.2695 ymin: 32.53432 xmax: -114.229 ymax: 41.99317
## Geodetic CRS:  WGS 84


Read and Write Rasters

To read in gridded data, use the rast() function. This will read in the Near Infrared (NIR) channel (B5)

b5 <- rast("./data/rs/LC08_044034_20170614_B5.tif")

b5
## class       : SpatRaster 
## dimensions  : 1245, 1497, 1  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 10N (EPSG:32610) 
## source      : LC08_044034_20170614_B5.tif 
## name        : LC08_044034_20170614_B5 
## min value   :            0.0008457669 
## max value   :            1.0124315023

When we print the rast object created, the second line of the output lists the dimensions of the data. Note that here, this has 1245 rows, 1497 columns and 1 layer. This also shows the resolution (30x30 m), the extent, the CRS and the a brief summary of the data.

We can write rast objects back to file using writeRaster() (I’ll bet you never thought it would be called that). You can write out to any format supported by GDAL. Here we write out to a TIFF format. You can see the full list of available formats for reading and writing by running the writeFormats() function in the console. We’ll use this again after having worked with the data.

writeRaster(b5, 
            filename = "./b5.tif",
            overwrite = TRUE)


Raster CRS

As with the sf objects, we can check the coordinate reference system of the file we just read. This does not print very well, but if you look at the end you’ll see a reference to the EPSG code (32610).

crs(b5)
## [1] "PROJCRS[\"WGS 84 / UTM zone 10N\",\n    BASEGEOGCRS[\"WGS 84\",\n        DATUM[\"World Geodetic System 1984\",\n            ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4326]],\n    CONVERSION[\"UTM zone 10N\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",-123,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",500000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"Engineering survey, topographic mapping.\"],\n        AREA[\"Between 126°W and 120°W, northern hemisphere between equator and 84°N, onshore and offshore. Canada - British Columbia (BC); Northwest Territories (NWT); Nunavut; Yukon. United States (USA) - Alaska (AK).\"],\n        BBOX[0,-126,84,-120]],\n    ID[\"EPSG\",32610]]"

If the CRS is not set, you can set it using the crs() function and an EPSG code. For example, the following code would set the CRS to WGS84 (don’t run this as the CRS is already defined)

crs(b5) <- "EPSG:4326"

Again, this should not be used to change the CRS, only set it. Note that crs is for rast objects, st_crs for vectors.


You can transform the CRS for a raster layer using project(). This can again use a EPSG type code.

b5_wgs84 <- project(b5, "EPSG:4326")

crs(b5_wgs84)
## [1] "GEOGCRS[\"WGS 84\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    USAGE[\n        SCOPE[\"Horizontal component of 3D system.\"],\n        AREA[\"World.\"],\n        BBOX[-90,-180,90,180]],\n    ID[\"EPSG\",4326]]"

We’ll keep the Landsat data in its UTM projection (10N), and reproject the CA place data to match:

ca_places <- st_transform(ca_places, 32610)


Basic Plotting

You can make a simple plot using the plot() function:

plot(b5, main = "Landsat 8 (B2)")

plot(st_geometry(ca_places), add = TRUE)


Summary Statistics

The function cellStats() can be used to calculate most summary statistics for a raster layer. So to get the mean global temperature (and standard deviation):

global(b5, mean)
##                              mean
## LC08_044034_20170614_B5 0.2693644
global(b5, sd)
##                               sd
## LC08_044034_20170614_B5 0.111683


Subset Rasters

If we want to use only a subset of the original raster layer, the function crop() will extract only the cells in a given region. This can be defined using another raster object or Spatial* object, or by defining an extent object:

# extent method
my_ext <- ext(c(xmin = 612500, 
                    xmax = 617500, 
                    ymin = 4196000,
                    ymax = 4201000))

b5_sub <- crop(b5, my_ext)

plot(b5_sub)

We can also use an sf object to crop the data. Here, we’ll extract the polygon for Bethel Island from the ca_places object, and use this to crop the raster:

bethel <- ca_places %>% 
  dplyr::filter(NAME == "Bethel Island")

b5_sub <- crop(b5, bethel)

plot(b5_sub)
plot(st_geometry(bethel), add = TRUE)

Note that crop subsets the original raster to the extent of Canada’s borders, rather than to the borders themselves. This is because rasters are always rectangular. You can ‘hide’ the values of raster cells outside of a polygon by using the mask function. The raster has to be rectangular, so this does not remove the cells outside the polygon. Rather, it sets their value to NA.

b5_sub <- mask(b5_sub, mask = bethel)

plot(b5_sub)
plot(st_geometry(bethel), add = TRUE)

Extract Data

Values can be extracted from individual locations (or sets of locations) using extract(). This can take a set of coordinates in matrix form, or use a Spatial* object. To get the reflectance value at 615000 E, 4199000 N:

extract(b5, cbind(615000,4199000))
##   LC08_044034_20170614_B5
## 1               0.2121801

You can also extract for multiple locations. Let’s generate a set of random points in Bethel Island, and then sample the reflectance value for each of these.

random_pnts <- st_sample(bethel, size = 20)

extract(b5, st_coordinates(random_pnts))
##    LC08_044034_20170614_B5
## 1               0.10249011
## 2               0.27481058
## 3               0.35472512
## 4               0.20094655
## 5               0.06612195
## 6               0.09268785
## 7               0.21324277
## 8               0.32954717
## 9               0.36669603
## 10              0.46638858
## 11              0.41353872
## 12              0.31954971
## 13              0.39768595
## 14              0.34171325
## 15              0.46441513
## 16              0.48172089
## 17              0.36231536
## 18              0.39916062
## 19              0.34785050
## 20              0.28463453

You can also extract values within a polygon, by replacing the point coordinates with a sf polygon:

b5_bethel <- extract(b5, bethel)

head(b5_bethel)
##   ID LC08_044034_20170614_B5
## 1  1              0.10251180
## 2  1              0.08084705
## 3  1              0.05189565
## 4  1              0.17219034
## 5  1              0.29502234
## 6  1              0.31579795

By default, this returns the value of all pixels within the polygon. By adding the fun= argument, you can easily calculate summary statistics:

extract(b5, bethel, fun = 'median')
##   ID LC08_044034_20170614_B5
## 1  1               0.3779947

Note that if the sf object has multiple polygons, it will return the summary statistic for each one. Let’s now extract the values for a different place (Oakley), and then we can compare them. We do this in a couple of steps. First get the polygon for Oakley, then extract using the two polygons combined.

oakley <- ca_places %>% 
  dplyr::filter(NAME == "Oakley")

b5_bethel_oakley <- extract(b5, rbind(bethel, oakley))

names(b5_bethel_oakley) <- c("ID", "B5")

We’ll visualize the difference with ggplot2, which shows higher NIR reflectance values for Bethel, likely indicating higher vegetation cover.

ggplot(b5_bethel_oakley, aes(x = B5, fill = as.factor(ID))) +
  geom_histogram(alpha = 0.7, position = 'identity')
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.


Raster Stacks

The b2 raster represents a single band from the Landsat 8 scene. More usefully, we can load several bands and combine them into a single object. Here, we’ll load the blue (B2), green (B3), red (B4), and infrared (B5) bands :

# Blue
b2 <- rast('data/rs/LC08_044034_20170614_B2.tif')
# Green
b3 <- rast('data/rs/LC08_044034_20170614_B3.tif')
# Red
b4 <- rast('data/rs/LC08_044034_20170614_B4.tif')
# Near Infrared (NIR)
b5 <- rast('data/rs/LC08_044034_20170614_B5.tif')

Now we’ll create a raster stack with all 4 of these:

s <- c(b5, b4, b3, b2)

s
## class       : SpatRaster 
## dimensions  : 1245, 1497, 4  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 10N (EPSG:32610) 
## sources     : LC08_044034_20170614_B5.tif  
##               LC08_044034_20170614_B4.tif  
##               LC08_044034_20170614_B3.tif  
##               LC08_044034_20170614_B2.tif  
## names       : LC08_04~0614_B5, LC08_04~0614_B4, LC08_04~0614_B3, LC08_04~0614_B2 
## min values  :    0.0008457669,      0.02084067,      0.04259216,       0.0748399 
## max values  :    1.0124315023,      0.78617686,      0.69246972,       0.7177562

The metadata now shows that this contains 4 layers. You can also make the stack directly by passing a list of file names:

filenames <- paste0('./data/rs/LC08_044034_20170614_B', 1:11, ".tif")

landsat <- rast(filenames)

landsat
## class       : SpatRaster 
## dimensions  : 1245, 1497, 11  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 10N (EPSG:32610) 
## sources     : LC08_044034_20170614_B1.tif  
##               LC08_044034_20170614_B2.tif  
##               LC08_044034_20170614_B3.tif  
##               ... and 8 more source(s)
## names       : LC08_~14_B1, LC08_~14_B2, LC08_~14_B3, LC08_~14_B4,  LC08_~14_B5,  LC08_~14_B6, ... 
## min values  :  0.09641791,   0.0748399,  0.04259216,  0.02084067, 0.0008457669, -0.007872183, ... 
## max values  :  0.73462820,   0.7177562,  0.69246972,  0.78617686, 1.0124315023,  1.043204546, ...

This now contains all bands representing reflection intensity in the following wavelengths: Ultra Blue, Blue, Green, Red, Near Infrared (NIR), Shortwave Infrared (SWIR) 1, Shortwave Infrared (SWIR) 2, Panchromatic, Cirrus, Thermal Infrared (TIRS) 1, Thermal Infrared (TIRS) 2.

If we now use the extract() function with this object, it will return values for all bands in the stack:

extract(landsat, cbind(615000,4199000))
##   LC08_044034_20170614_B1 LC08_044034_20170614_B2 LC08_044034_20170614_B3
## 1               0.1691325               0.1547111               0.1395956
##   LC08_044034_20170614_B4 LC08_044034_20170614_B5 LC08_044034_20170614_B6
## 1               0.1370149               0.2121801               0.1594604
##   LC08_044034_20170614_B7 LC08_044034_20170614_B8 LC08_044034_20170614_B9
## 1               0.1526508              0.09177701             0.001236123
##   LC08_044034_20170614_B10 LC08_044034_20170614_B11
## 1                 311.6736                  308.671
par(mfrow = c(2,2))
plot(b2, main = "Blue", col = gray(0:100 / 100))
plot(b3, main = "Green", col = gray(0:100 / 100))
plot(b4, main = "Red", col = gray(0:100 / 100))
plot(b5, main = "NIR", col = gray(0:100 / 100))

The values in each layer range from 0 to 1, and the same scale has been used for each band, showing clearly the difference in reflectance for this set of wavelengths. For example, vegetation reflects more energy in NIR than other wavelengths and thus appears brighter. In contrast, water absorbs most of the energy in the NIR wavelength and it appears dark

Composite images

The bands can be combined to form composite images. Here, we’ll use the red, green and blue bands make a true color image. This uses the concatenation function (c()), and the order of the bands is important (R, then G, then B):

landsatRGB <- c(b4, b3, b2)

plotRGB(landsatRGB, stretch = "lin")

We can also make a false color composite with the NIR, red and green bands, where bright reds indicate vegetation cover:

Another popular image visualization method in remote sensing is known “false color” image in which NIR, red, and green bands are combined. This representation is popular as it makes it easy to see the vegetation (in red).

landsatFCC <- c(b5, b4, b3)
plotRGB(landsatFCC, stretch="lin")

# Raster algebra

terra makes it easy to carry out simple raster algebraic operations. Eahc band or layer is treated a 2D array which makes it possible to add, subtract, multiply, divide, etc. As a simple example here, we can calculate the NDVI for the Landsat scene as

\[ NDVI = (NIR - R) / (NIR + R) \]

NIR is band 5 and red is band 4:

ndvi <- (b5 - b4) / (b5 + b4)

plot(ndvi, col=rev(terrain.colors(10)), main = "NDVI")

We can also make a quick histogram of values to look for any outliers:

hist(ndvi, main = "NDVI")
## Warning: [hist] a sample of54% of the cells was used

We can then use the results to carry out some simple classification.

  • Vegetation (NDVI > 0.4). The clamp function masks all value otuside of a range (here below 0.4)
veg = clamp(ndvi, lower=0.4, values=FALSE)

plot(veg)

  • Croplands (corresponding to the peak in NDVI values):
crops = ndvi > 0.25 & ndvi < 0.3 

plot(crops)

And water bodies (NDVI < 0):

water = ndvi < 0

plot(water)

Mapping packages

In the previous lab, we looked at making maps in base R (with the plot() function and an sf object) as well as using ggplot2. More recently, there has been something of an explosion in mapping packages, including interactive maps. We’ll look here at two of these: tmap and leaflet. Other useful ones include mapsf highcharter, mapview and mapdeck.

library(sf)

## Census tracts for Salt Lake County with population density
tracts <- st_read("./data/slc_tract/slc_tract_2015.shp")
## Reading layer `slc_tract_2015' from data source 
##   `/Users/u0784726/Dropbox/Data/devtools/ugic-intro-r/data/slc_tract/slc_tract_2015.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 212 features and 17 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -112.2602 ymin: 40.41417 xmax: -111.5532 ymax: 40.92185
## Geodetic CRS:  NAD83
## Salt Lake light rail tracks
lightrail <- st_read("./data/LightRail_UTA/LightRail_UTA.shp")
## Reading layer `LightRail_UTA' from data source 
##   `/Users/u0784726/Dropbox/Data/devtools/ugic-intro-r/data/LightRail_UTA/LightRail_UTA.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 38 features and 5 fields
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: 413241.3 ymin: 4486347 xmax: 429455.1 ymax: 4515254
## Projected CRS: NAD83 / UTM zone 12N
## Salt Lake light rail stations
stations <- st_read("./data/LightRailStations_UTA/LightRailStations_UTA.shp")
## Reading layer `LightRailStations_UTA' from data source 
##   `/Users/u0784726/Dropbox/Data/devtools/ugic-intro-r/data/LightRailStations_UTA/LightRailStations_UTA.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 57 features and 6 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 413245.8 ymin: 4486448 xmax: 429389.5 ymax: 4515197
## Projected CRS: NAD83 / UTM zone 12N
## Landsat images for California
# Blue
b2 <- rast('data/rs/LC08_044034_20170614_B2.tif')
# Green
b3 <- rast('data/rs/LC08_044034_20170614_B3.tif')
# Red
b4 <- rast('data/rs/LC08_044034_20170614_B4.tif')
## Places for California
ca_places <- st_read("./data/ca_places/ca_places.shp")
## Reading layer `ca_places' from data source 
##   `/Users/u0784726/Dropbox/Data/devtools/ugic-intro-r/data/ca_places/ca_places.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1618 features and 16 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -124.2695 ymin: 32.53432 xmax: -114.229 ymax: 41.99317
## Geodetic CRS:  WGS 84
## W North America Climate
wna_climate <- read.csv("./data/WNAclimate.csv")
## Natural Earth country shapefiles 
countries <- st_read("./data/ne_50m_admin_0_countries/ne_50m_admin_0_countries.shp")
## Reading layer `ne_50m_admin_0_countries' from data source 
##   `/Users/u0784726/Dropbox/Data/devtools/ugic-intro-r/data/ne_50m_admin_0_countries/ne_50m_admin_0_countries.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 241 features and 94 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -180 ymin: -89.99893 xmax: 180 ymax: 83.59961
## Geodetic CRS:  WGS 84
## Natural Earth country rivers

tmap

tmap makes thematic maps. It works in a very similar way to ggplot2, by building a series of layers and map geometries and elements. In general, we start by using tm_shape() to identify the spatial object to be used, and add geometries are added, including filled polygons, borders and symbols. Finally, we can add legends, scale bars, etc. Note that you will need to install this package if you haven’t done so already (install.packages("tmap")).

Example 1: Salt Lake County

We’ll first make a map with a mixture of vector layers, including census tracts, light rail tracks and stations. Start by loading these

library(sf)

## Census tracts for Salt Lake County with population density
tracts <- st_read("./data/slc_tract/slc_tract_2015.shp")
## Reading layer `slc_tract_2015' from data source 
##   `/Users/u0784726/Dropbox/Data/devtools/ugic-intro-r/data/slc_tract/slc_tract_2015.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 212 features and 17 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -112.2602 ymin: 40.41417 xmax: -111.5532 ymax: 40.92185
## Geodetic CRS:  NAD83
## Salt Lake light rail tracks
lightrail <- st_read("./data/LightRail_UTA/LightRail_UTA.shp")
## Reading layer `LightRail_UTA' from data source 
##   `/Users/u0784726/Dropbox/Data/devtools/ugic-intro-r/data/LightRail_UTA/LightRail_UTA.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 38 features and 5 fields
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: 413241.3 ymin: 4486347 xmax: 429455.1 ymax: 4515254
## Projected CRS: NAD83 / UTM zone 12N
## Salt Lake light rail stations
stations <- st_read("./data/LightRailStations_UTA/LightRailStations_UTA.shp")
## Reading layer `LightRailStations_UTA' from data source 
##   `/Users/u0784726/Dropbox/Data/devtools/ugic-intro-r/data/LightRailStations_UTA/LightRailStations_UTA.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 57 features and 6 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 413245.8 ymin: 4486448 xmax: 429389.5 ymax: 4515197
## Projected CRS: NAD83 / UTM zone 12N

If we check the CRS, you’ll see that tracts has a different one, so will need to be reprojected:

st_crs(tracts)$epsg
## [1] 4269
st_crs(lightrail)$epsg
## [1] 26912
tracts <- st_transform(tracts, st_crs(lightrail))

Now let’s start mapping things out with tmap. First, let’s make a simple map showing the polygon outlines using tm_borders(). This is the standard tmap method: first define the object that you want to plot with tm_shape, then define what you want to plot.

library(tmap)

tm_shape(tracts) + 
  tm_borders()

The function tm_fill() will then fill these using one of the variables in the tract data set. We’ll use the population density (density). Note that this automatically adds a legend within the frame of the figure:

tm_shape(tracts) + 
  tm_fill("density") +
  tm_borders()

The color scale can be changed by setting the palette argument in tm_fill(). This includes the ColorBrewer scales described in the previous module, and the algorithm used to define the intervals. For example, to use the ‘Greens’ palette with percentile breaks. We also change the color of the tract boundaries to lightgray:

tm_shape(tracts) + 
  tm_fill("density", palette = "Greens", style = "quantile") +
  tm_borders("lightgray")

You can reverse the color scale by prepending a - (-Greens). If you want to see the full set of palettes that you can use, install the tmaptools package and run the following code (you may also need to install shinyjs):

tmaptools::palette_explorer()

Let’s now add another layer. We’ll overlay the light rail track on this map. As this is a different object in R (lightrail), we first need to use tm_shape to indicate that we are using it, then add tm_lines() to show the routes. We’ll make the lines a little thicker with lwd and change the type to dashed.

tm_shape(tracts) + 
  tm_fill("density", palette = "Greens", style = "quantile") +
  tm_borders("lightgray") +
  tm_shape(lightrail) +
  tm_lines(lwd = 2, lty = 'dashed', col = "darkorange")

Note that if you use a varible name for the color, it will thematically map the lines (I’ve removed the polygon fill to make this clear):

tm_shape(tracts) + 
  #tm_fill("density", palette = "Greens", style = "quantile") +
  tm_borders("lightgray") +
  tm_shape(lightrail) +
  tm_lines(lwd = 4, col = "ROUTE")

We can now add the stations (using tm_shape() to select the correct object to plot):

tm_shape(tracts) + 
  tm_fill("density", palette = "Greens", style = "quantile") +
  tm_borders("lightgray") +
  tm_shape(lightrail) +
  tm_lines(lwd = 2) +
  tm_shape(stations) +
  tm_dots(size = 0.25, shape = 23)

Now let’s add some details to the map. We’ll add a graticule before the polygons are plotted, a compass and a scalebar. We’ll also use the tm_layout function to add a title and move the legend to a new position. (This function has a lot of options for changing your map - it’s well worth checking the help page.)

tm_shape(tracts) + 
  tm_graticules(col = "lightgray") + 
  tm_fill("density", title = "Popn density", palette = "Greens", style = "quantile") +
  tm_borders("lightgray") +
  tm_shape(lightrail) +
  tm_lines(lwd = 2) +
  tm_shape(stations) +
  tm_dots(size = 0.25, shape = 23) +
  tm_compass(position = c("left", "bottom")) +
  tm_scale_bar(position = c("right", "top")) +
  tm_layout(main.title = "Salt Lake County Light Rail", 
            legend.outside = TRUE,
            legend.outside.position = c("left"))

The map can be clipped to smaller regions by using the bbox argument in the first tm_shape() function. In this code, we’ll extract a set of census tracts that correspond to downtown Salt Lake. We’ll then use the the bounding box of these to crop the map (I’ve removed the fill). We’ll also add the station name using tm_text():

tracts_sub = tracts %>%
  dplyr::filter(TRACTCE %in% c(102500, 102600, 114000))

tm_shape(tracts, bbox = st_bbox(tracts_sub)) + 
  tm_borders("lightgray") +
  tm_shape(lightrail) +
  tm_lines(lwd = 2) +
  tm_shape(stations) +
  tm_dots(size = 0.25, shape = 23) +
  tm_text("STATIONNAM", ymod = -1, bg.color = "white", size = 0.8)

Finally, tmap allows you to make interactive maps fairly simply with the tmap_mode() function. Setting this to view will make an interactive map, to plot will make a static map (like the ones we have made so far). Here we’ll remake the full map as an interactive map:

## Set interactive
tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(tracts) + 
  tm_fill("density", title = "Popn density", palette = "Greens", style = "quantile", id = "TRACTCE") +
  tm_borders("lightgray") +
  tm_shape(lightrail) +
  tm_lines(lwd = 2) +
  tm_shape(stations) +
  tm_dots(size = 0.25, shape = 23)
## Symbol shapes other than circles or icons are not supported in view mode.

We’ll reset back to static maps for the rest of this exercise:

tmap_mode("plot")
## tmap mode set to plotting

Example 2: W N America climate

In this example, we’ll map out the climate data from western North America. We’ll go a bit faster in this example and not explain every step. First read in the data and create an sf object:

sf_use_s2(FALSE)
## Spherical geometry (s2) switched off
wna_climate <- read.csv("./data/WNAclimate.csv")

wna_climate <- st_as_sf(wna_climate, 
                        coords = c("LONDD", "LATDD"),
                        crs = 4326)

Make a quick thematic map of average july temperature:

tm_shape(wna_climate) +
  tm_symbols(col = "Jul_Tmp")

Next, we’ll download country outlines and river centerlines from Natural Earth using the rnaturalearth package (you will need to install this). This will download layers from the Natural Earth website - there are a variety of these available, including country and state boundaries, physical features (rivers, lakes, etc) and cultural features (place names, roads, etc). There are three levels of resolution: 1:10m, 1:50m and 1:110m. To download a layer, you need to specify the resolution (scale) and the name of the layer. The full set can be found here

library(rnaturalearth)

countries50 <- ne_download(scale = 50, type = "admin_0_countries")
## Reading layer `ne_50m_admin_0_countries' from data source 
##   `/private/var/folders/ql/nw995vq50pq3dlrxhk7mm4_40000gq/T/RtmpiMUnCe/ne_50m_admin_0_countries.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 242 features and 168 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -180 ymin: -89.99893 xmax: 180 ymax: 83.59961
## Geodetic CRS:  WGS 84
rivers50 <- ne_download(scale = 50, type = "rivers_lake_centerlines", category = "physical")
## Reading layer `ne_50m_rivers_lake_centerlines' from data source 
##   `/private/var/folders/ql/nw995vq50pq3dlrxhk7mm4_40000gq/T/RtmpiMUnCe/ne_50m_rivers_lake_centerlines.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 478 features and 36 fields (with 1 geometry empty)
## Geometry type: MULTILINESTRING
## Dimension:     XY
## Bounding box:  xmin: -165.2439 ymin: -50.24014 xmax: 176.3258 ymax: 73.3349
## Geodetic CRS:  WGS 84

Now we can put this all together as a series of layers. We’ll make the color palette continuous and use a red-blue palette.

tm_shape(countries50, bbox = st_bbox(wna_climate)) +
  tm_borders() +
  tm_shape(rivers50) +
  tm_lines("lightblue", lwd = 2) +
  tm_shape(wna_climate) +
  tm_symbols(col = "Jul_Tmp", palette = "-RdBu", alpha = 0.75,
             style = "cont", title.col = "degC")
## Warning: The shape rivers50 contains empty units.

Finally, we’ll make two maps (july temperature and annual precipitation), and use tmap_arrange() to make a two panel plot:

m1 = tm_shape(countries50, bbox = st_bbox(wna_climate)) +
  tm_borders() +
  tm_shape(rivers50) +
  tm_lines("lightblue", lwd = 2) +
  tm_shape(wna_climate) +
  tm_symbols(col = "Jul_Tmp", palette = "-RdBu", alpha = 0.75,
             style = "cont", title.col = "degC")

m2 = tm_shape(countries50, bbox = st_bbox(wna_climate)) +
  tm_borders() +
  tm_shape(rivers50) +
  tm_lines("lightblue", lwd = 2) +
  tm_shape(wna_climate) +
  tm_symbols(col = "annp", palette = "BuPu", alpha = 0.75,
             style = "cont", title.col = "mm/yr")

tmap_arrange(m1, m2)
## Warning: The shape rivers50 contains empty units.

## Warning: The shape rivers50 contains empty units.

## Warning: The shape rivers50 contains empty units.

## Warning: The shape rivers50 contains empty units.

Example 3: Raster images

In the final example, we’ll use the Landsat data with tmap. The main function to map raster data is tm_raster(). We’ll use it here to show the ndvi layer made earlier.

tm_shape(ndvi) +
  tm_raster()
## stars object downsampled to 1096 by 912 cells. See tm_shape manual (argument raster.downsample)
## Variable(s) "NA" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.

We’ll update this a little by clamping the NDVI values at a lower limit of 0, changing the color palette and making it continuous. We’ll also add the California places added as an overlay.

tm_shape(clamp(ndvi, lower = 0)) +
  tm_raster(palette = "Greens", style = "cont", title = "NDVI") +
  tm_shape(ca_places) +
  tm_borders(lwd = 2) +
  tm_layout(main.title = "Landsat 8 (2017/06/14)", 
            legend.position = c("left", "top"), 
            legend.bg.color = "white", legend.bg.alpha = 0.7)
## stars object downsampled to 1096 by 912 cells. See tm_shape manual (argument raster.downsample)

A second function (tm_rgb()) allows the creation of color composites. We’ll use the red, green and blue bands to make the true color composite. The maximum reflectance value in these files is 1.0, so we need to set this for scaling.

tm_shape(landsat) + 
  tm_rgb(r = 4, g = 3, b = 2, max.value = 1)
## stars object downsampled to 1096 by 912 cells. See tm_shape manual (argument raster.downsample)
## Warning: Raster values found that are outside the range [0, 1]

When we made this figure earlier, it was a lot brighter. This is because the previous function ‘stretched’ the range of values for each layer. We can mimic that here by first using the stretch() function from terra, then repeating the plot:

landsat_stretch = stretch(landsat, minv = 0, maxv = 1, 
                          minq = 0.01, maxq = 0.99)

tm_shape(landsat_stretch) + 
  tm_rgb(r = 4, g = 3, b = 2, max.value = 1)
## stars object downsampled to 1096 by 912 cells. See tm_shape manual (argument raster.downsample)

Leaflet

The tmap mode function allows you to quickly make an interacitve map. The resulting map is based on the well-known Leaflet javascript library, which offers a much greater range of flexibility in making these maps. Fortunately, there is an R library (leaflet) that allows access to the Leaflet API allowing you to build quite complex maps from R. We’ll remake the map of Salt Lake County and the light rail routes using this now. First install the library (install.packages("leaflet")), and load it (and a helper library for working with html formats):

library(leaflet)
library(htmltools)

leaflet operates in a similar way to ggplot2 and tmap in building a map by layers. Layers are added using the pipe operator that we encountered in an earlier lab (%>%). There are three basic ingredients that are required:

  • A blank map canvas from the leaflet() function
  • A set of background tiles
  • One or more markers or polygons to display

Here’s a really simple example, where we use the default tiles from OpenStreetMap and drop a marker on the University of Utah:

m <- leaflet() %>%
  addTiles() %>%  # Add default OpenStreetMap map tiles
  addMarkers(lng=-111.8421, lat=40.7649, popup="The University of Utah")
m  # Print the map

You can also try changing the Provider - the source of the background map. For example, this will use ESRI map (the map is not shown here but will appear in RStudio when you run this):

m <- leaflet() %>%
  addProviderTiles("Esri.WorldStreetMap") %>%  
  addMarkers(lng=-111.8421, lat=40.7649, popup="The University of Utah")
m

ESRI imagery:

m <- leaflet() %>%
  addProviderTiles("Esri.WorldImagery") %>%  # Add default OpenStreetMap map tiles
  addMarkers(lng=-111.8421, lat=40.7649, popup="The University of Utah")
m

You can find a full set of these providers here: https://leaflet-extras.github.io/leaflet-providers/preview/

We’ll now add the census tracts and light rail sf objects. First, these need to be reprojected to WGS84 long-lat (EPSG code 4326) to work with Leaflet.

tracts_ll <- st_transform(tracts, crs = 4326)
stations_ll <- st_transform(stations, crs = 4326)
lightrail_ll <- st_transform(lightrail, crs = 4326)

Now we can add these to a dark basemap. First the tracts with addPolygons, then the routes (addPolylines) and stations (addMarkers). Note that the stations include a label taken from the station name column in the sf object. This wil appear when you hover your cursor over a marker.

leaflet() %>%
  # add a dark basemap
  addProviderTiles("CartoDB.DarkMatter") %>%
  # add the polygons of the clusters
  addPolygons(
    data = tracts_ll,
    color = "#E2E2E2",
    opacity = 1, # set the opacity of the outline
    weight = 1, # set the stroke width in pixels
    fillOpacity = 0.2 # set the fill opacity
  ) %>%
  addPolylines(
    data = lightrail_ll
  ) %>%
  addMarkers(
    data = stations_ll,
    label = ~htmlEscape(STATIONNAM)
  )

Finally, we’ll add some extra details. First, we’ll set colors for each station by the light rail line it lies on. This is a little complicated - below we:

  • Get the list of unique station names
  • Make up a color palette from the RColorBrewer package (you may need to install this)
  • Convert the station names to an integer code (so the first name alphabetically will be 1, etc)
  • Finally use the integer code to assign the matching color back to the sf object
unique(stations_ll$LINENAME)
## [1] "University"            "N/S TRAX"              "Med Center"           
## [4] "Intermodal Hub"        "Airport"               "Draper"               
## [7] "Mid-Jordan"            "West Valley"           "Sugar House Streetcar"
line_pal <- RColorBrewer::brewer.pal(9, "Set1")
line_no <- as.numeric(as.factor(stations_ll$LINENAME))

stations_ll$LINECOL <- line_pal[line_no]

Next we’ll make a popup window for each station. This will appear when the marker is clicked (as opposed to the label which appears when a cursor hovers over the marker). This has to be formatted using html. Here, we use R’s paste() function to stick together html tags with the station name, address, etc. The <b> tags make the labels bold, and the <br> tags add a line break. You can do a lot more here, for example, including images or URLs in each popup.

station_popup = paste0(
  "<b>Station: </b>",
  stations_ll$STATIONNAM,
  "<br>",
  "<b>Line Name: </b>",
  stations_ll$LINENAME,
  "<br>",
  "<b>Park n Ride: </b>",
  stations_ll$PARKNRIDE,
  "<br>",
  "<b>Address: </b>",
  stations_ll$ADDRESS
  
)

With all this in place, we can rebuild our map. The main changes here are that we use a circle marker as this has an option to change color, and include the popup information in the marker.

leaflet() %>%
  # add a dark basemap
  addProviderTiles("Esri.WorldStreetMap") %>%
  # add the polygons of the clusters
  addPolygons(
    data = tracts_ll,
    color = "#E2E2E2",
    # set the opacity of the outline
    opacity = 1,
    # set the stroke width in pixels
    weight = 1,
    # set the fill opacity
    fillOpacity = 0.2
  ) %>%
  addPolylines(
    data = lightrail_ll
  ) %>%
  addCircleMarkers(
    data = stations_ll,
    color = ~LINECOL,
    label = ~htmlEscape(STATIONNAM),
    popup = station_popup
  )